BACKGROUND OF THE INVENTION 

Field of the Invention 

The present invention relates to a method allowing determination of a velocity 
model of seismic waves picked up by receivers coupled with an underground 
5 formation, from multi-offset records of these waves. 

DESCRIPTION OF THE PRIOR ART 

Seismic reflection surveys are widely used in petroleum exploration, notably to 
produce images of the subsurface from the information contained in the waves 
propagated and reflected from the geologic discontinuities of the subsoil. 

10 More precisely, imagery methods use the kinematic information associated 

with the major seismic reflections (that is the travel times of the waves reflected on 
the main discontinuities of the subsoil) to determine a velocity macromodel of the 
subsoil, which will be used to convert the temporal seismic records to a depth image 
of the subsurface. 

15 Access to the kinematic information necessary for determination of the 

velocity model requires interpretation of the seismic events in the seismic multi-offset 
records. Now, seismic multi-offset data is generally characterized by a bad signal-to- 
noise ratio, hence the failure of conventional automatic picking and the necessity of 
carrying out a long and costly manual picking of the seismic data. Besides, for 3D 

20 seismic surveys (currently predominant in relation to 2D surveys), the quantity of 
data to be interpreted is a 4D volume that can reach sizes of the order of one 
hundred gigabytes, or even of the order of a terabyte, which makes the interpretative 
task even longer and tedious. 



In order to avoid this stage of seismic multi-offset collection interpretation, 
geophysicists have developed methods based on approximations of the geometry of 
the seismic events in the multi-offset collections. To establish these approximations, 
these methods put forward hypotheses on the subsoil complexity, hypotheses which 
5 can relate to the geometry of the geologic discontinuities of the subsoil and to the 
layer velocity variations. The method described by Taner and Kohler (1969) can for 
example be mentioned, which assumes a stratified medium of horizontal plane 
layers with homogeneous and isotropic layer velocities, as well as the method 
described by Levin (1971) which extends the previous method to sloping layers. 
10 Other variants have been proposed, but generally the existing methods are based on 
a hyperbolic hypothesis of the travel time curve in the multi-offset collections. This 
hypothesis is however violated from the moment that the geologic discontinuities of 
the subsoil are no longer plane and/or that the layer velocities exhibit lateral 
variations. 

15 SUMMARY OF THE INVENTION 

The invention is a method for best approaching the travel times of the multi- 
offset collections in cases where the geologic discontinuities of the subsoil have any 
geometry and where the layer velocities are moderately variable laterally, while 
limiting the human time required for interpretation of the seismic records. The travel 
20 times obtained are then treated by means of a prestack kinematic inversion method 
taking account of complex kinematics, such as travel time tomography, described for 
example in the following document: 

Bishop, T. et al., 1985, "Tomographic Determination of Velocity and Depth in 
Laterally Varying Media", Geophysics, 50 No. 6, 903-923. 



The method can also apply in cases where the layer velocities are laterally 
variable, but iteratively, the velocity model after updating at a given iteration by 
applying a prestack kinematic inversion method such as the aforementioned 
reflection tomographic method serving as the initial model for the next iteration. 

5 The method of the invention allows easier access to the prestack kinematic 

information associated with the events contained in the seismic records, which is 
combined with a prestack kinematic inversion technique for imagery of the geologic 
interfaces of the subsoil. 

What is referred to as "prestack kinematic information" is the travel times 
10 associated with the reflections recorded by source-pickup pairs located at variable 
distances from one another. Contrasting of "easier access" with manual picking of 
the seismic events in the multi-offset collections is performed. A "prestack kinematic 
inversion method" is understood to be a method allowing, from the kinematic 
information extracted from the seismic records (and not from an approximation of 
15 this kinematic information), to find the geometry of the associated reflectors and the 
layer velocities. 

The method according to the invention allows determination a velocity model 
of seismic waves picked up by seismic receivers coupled with an underground 
formation, in response to the emission of seismic waves in the subsoil by a seismic 
20 source, after reflection by geologic interfaces of the formation, from multi-offset 
records of these waves. 

For each seismic event located in the records and for each layer delimited by 
interfaces, the method comprises at least the following steps: 
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a) from the prestack seismic records, constructing an iso-offset collection 
from which kinematic information or travel times associated with the event are 
extracted; 

b) selecting a velocity range around a reference velocity in the layer, that is 
5 sampled with a predetermined interval; 

c) for each velocity sample, applying an inversion technique at a fixed 
velocity so as to determine, from the travel times extracted from the iso-offset 
collection, a geometry of the interface for the velocity sample concerned in order to 
obtain a series of interface/velocity pairs for the event; 

10 d) calculating kinematic information associated with each interface/velocity 

pair obtained, for source-receiver pairs corresponding to multi-offset collections 
existing in the seismic records; 

e) for each interface/velocity pair and for each multi-offset collection selected, 
evaluating the coherence between the multi-offset travel times thus calculated and 

15 the seismic records, and selecting for each multi-offset collection the travel time 
curve showing maximum coherence with the seismic records; 

f) applying a prestack kinematic inversion method using the multi-offset travel 
times obtained for all the multi-offset collections selected, in order to determine the 
geometry and the velocity of the layer considered; and 

20 g) iterating n times (n ^ 0) steps a) to f) by considering on each iteration the 

velocity model obtained during the previous iteration as the reference model to 
define the reference velocity of the new iteration. 



5 



According to an implementation mode suitable in cases where the velocity 
range selected is not precise enough at the end of either the previous iteration or of 
step a) to f), step g) is carried out n times with n ^ 1 . 

According to another implementation mode suitable in cases where the 
5 velocity distribution varies greatly laterally and/or in cases where no sufficiently 
precise a priori knowledge of the velocity distribution in the layer considered is 
available, step g) is carried out on offset ranges that are increasingly greater as 
iterations progress. 

According to another implementation mode suitable in cases where the 
10 velocity distribution varies greatly laterally and/or in cases where no sufficiently 
precise a priori knowledge of the velocity distribution in the layer considered is 
available, step g) is carried out on multi-offset collection grids that are increasingly 
finer as iterations progress. 

According to an implementation mode suitable in cases where the interface 
15 geometry which is sought generates triplications, step e) is carried out by 

considering the triplication branches in the multi-offset collections independently of 
one another. For example, ray tracing and inversion tools allowing taking into 
account of the multi-valuated arrivals are used. 

According to an implementation example, a zero-offset or a near-offset iso- 
20 offset collection is preferably constructed in step a). 

According to another implementation example, a fixed-velocity kinematic 
inversion technique such as a map migration is advantageously applied in step c). 



6 



According to another implementation example, the kinematic information is 
calculated in step d) by tracing multi-offset rays on the interface of each interface- 
velocity pair. 

According to another implementation example, step f) is advantageously 
5 carried out by applying a prestack kinematic inversion method such as a prestack 
travel time tomography. 

In its most general definition corresponding to the case where the number n of 
iterations to be carried out may possibly be zero, the method applies for example if 
the velocity distribution in each layer is hardly variable laterally and/or in the case 
10 where a sufficiently precise a priori on the velocity distribution in the layer considered 
is available. 

In a more restrictive sense where the number n of iterations to be carried out 
is at least 1 , the method applies for example in the case where the velocity 
distribution in the layers of the geologic formation is greatly variable laterally and/or 
15 in the case where no sufficient a priori on the velocity distribution in the layer 
considered is available. 

BRIEF DESCRIPTION OF THE DRAWINGS 

Other features and advantages of the method according to the invention will 
be clear from reading the description hereafter of non limitative embodiment 
20 examples, with reference to the accompanying drawings wherein: 

Figures 1a to 1c respectively show an application of the method to a common 
receiver type multi-offset collection (Fig. 1a), to a common receiver collection on 
which are superposed the multi-offset travel times obtained with various layer 



velocity values (Fig. 1b) and to a common receiver collection on which is superposed 
the predicted multi-offset travel time curve showing the greatest coherence with the 
seismic event considered (Fig. 1c). 

DETAILED DESCRIPTION OF THE INVENTION 

5 Seismic records are obtained by means of a seismic device comprising a 

seismic source emitting seismic waves propagated in the subsoil, a set of seismic 
receivers are coupled with the medium which pick up the waves reflected by the 
subsoil discontinuities in response to the waves emitted, and a laboratory is used for 
recording the picked up seismic signals. 

10 I - Standard case 

A first implementation mode of the method is described, in the case where the 
velocity distribution in the geologic formation preceding the layer being examined is 
known or has been estimated otherwise, the geometry of the layer in question 
generating no triplications, the velocity to be determined being slightly variable 
15 laterally or a sufficiently precise a priori on this velocity being available. The 
following steps are carried out: 

1-1 From the prestack or multi-offset seismic records, a collection of constant 
offset, is constructed preferably zero-offset (in reality an approximation of a zero- 
offset collection, known as a stack, characterized by a better signal-to-noise ratio 
20 than the individual iso-offset sections) or near-offset (corresponding to the lowest 
offset of the multi-offset records). 

If a coherent and sufficiently complete interpretation of the seismic event in 
progress can be carried out in this collection, picking is directly performed in this 



collection. In the opposite case, in order to gain access to the travel times 
associated with the event in question, a detour can be made via the time or depth 
migrated domain, followed by a time or depth demigration, as in the method 
described in: 

5 Ehinger, A., and Lailly, P., 1995, Velocity Model Determination by the SMART 
Method, Part 1: Theory: 65 th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded 
Abstracts, pp. 739-742. 

I-2 A reference velocity is selected for the layer in question, as well as an 
uncertainty of this reference velocity and a sampling interval for the velocity range 

10 thus formed. Thus a series of velocity samples for the layer in question is obtained; 
the reference velocity can be any velocity (that is, laterally and vertically variable); in 
particular, if a velocity distribution obtained otherwise is available, this velocity can 
be used as the reference velocity; if no a priori information is available for the 
velocity of the layer, a laterally and vertically homogeneous reference velocity can, 

15 for example, be taken. 

1-3 For each velocity sample, a fixed-velocity type kinematic inversion 
technique such as a map migration for example is applied in order to determine, from 
the travel times extracted from the constant-offset section, the geometry of the 
interface for the velocity sample in question. Thus a series of interface/velocity pairs 
20 for the event in question is obtained. A description of map migration can for example 
be found in the following reference: 

Yilmaz, O., 2001 , Seismic Data Analysis - Processing, Inversion and Interpretation of 
Seismic Data: Society of Exploration Geophysicists. 
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It may be impossible to find an interface model explaining the travel times for 
the fixed velocity in progress (impossibility of reaching convergence). For a higher 
efficiency of the method, it is possible to select only interface/velocity pairs for which 
the difference between the travel times observed and the travel times calculated on 
5 the model obtained after convergence is below a certain threshold. 

It can be noted that this fixed-velocity inversion can also be carried out by 
means of any other equivalent kinematic technique. 

1-4 The kinematic information associated with each interface/velocity pair 
obtained for example by means of multi-offset ray tracing on the interface of each 
10 interface/velocity pair thus obtained and for source-receiver pairs corresponding to 
multi-offset collections existing in the seismic records is calculated. Common 
midpoint collections are preferably used, in particular if the real velocity shows lateral 
variations that are not taken into account in the reference model. 

1-5 For each interface/velocity pair and for each multi-offset collection 
15 selected, the coherence between the multi-offset travel times thus calculated and the 
seismic records is evaluated. A technique referred to as semblance calculation 
technique can for example be used, wherein the energy encountered in the seismic 
data along a curve that corresponds here to the multi-offset travel times calculated 
by ray tracing is summed. In an ideal case, there will be a velocity sample for which 
20 the predicted travel time curve is perfectly superposed on the seismic event 

considered, which will translate into a maximum coherence measurement. It can be 
noted that this evaluation can be performed on a limited offset range, defined by the 
user or automatically, for example by examining the curve showing the evolution of 
the coherence as a function of the offset range. 



For each multi-offset collection, the travel time curve which best matches the 
seismic records is selected, and the multi-offset travel times forming this curve are 
reserved. It can be noted that the travel time curve exhibiting the best coherence 
can be adjusted if necessary by seeking in a predetermined vertical window the 
5 closest amplitude maximum (or minimum, depending on the polarity of the event 
considered) for each seismic trace of the multi-offset collection considered. 

1-6 The multi-offset travel times thus collected for all the multi-offset 
collections selected are then injected into a prestack kinematic inversion method 
such as a travel time tomography in order to determine the geometry and the velocity 
10 of the layer in question. 

I- 7 The previous stages may have to be iterated if the coherence with the 
seismic records is not considered to be globally satisfactory by the operator. 

II - Particular cases 

II- 1 In cases where the velocity distribution in the layers of the geologic 

15 formation is greatly variable laterally and/or where no sufficiently precise a priori on 
the velocity distribution in the layer considered is available, the method according to 
the invention comprises the following steps: 

The previous stages are applied iteratively, the travel times collected during 
an iteration being used for updating the velocity model by means of a prestack 
20 kinematic inversion method such as travel time tomography, a model which is then 
used as the input velocity model for a new iteration of the method. 

It can be noted that, in order to be free from the effects of the lateral velocity 
variations that might have a smaller wavelength than the maximum offset of the 



multi-offset collections, the offset range considered by the method can be widened 
as the iterations progress. More precisely, the method is initiated on a limited offset 
range, then the velocity distribution found from this limited offset range is used as the 
reference model for a new iteration of the method during which a wider offset range 
5 will be considered, and so forth. It can be noted that the offset range considered at a 
given iteration can be selected for example by examining the curve showing the 
evolution of the coherence as a function of the offset range. 

For the cases where the velocity distribution varies greatly laterally, it is 
possible to adopt furthermore a multi-grid approach, more precisely, the method is 
10 initiated on a loose grid of multi-offset collections (in order to determine the great 
wavelengths of the velocity variations), then multi-offset collection grids that are 
increasingly finer as iterations progress are taken into account. 

11-2 In the case where the sought interface geometry generates triplications 
and if it is also desired to determine the zones of the layer considered generating 
15 these triplications, it is first and foremost necessary to have inversion methods such 
as (fixed-velocity and prestack) kinematic inversion methods for example, allowing to 
take account of multi-valuated arrivals. Such a method is described in the following 
reference: 

Delprat-Jannaud, F. and Lailly, P., (1995), How to Handle Multiple Arrivals Journal of 
20 Geophysical Research, 100, No. B1, 703-715. 

The input data of this type of method are the multi-valuated travel times 
associated with the seismic events considered, and the ray parameters associated 
with these travel times. In order to obtain the multi-valuated travel times and the ray 
parameters associated with the constant-offset data, a migration-demigration loop 



can be used (as described in the aforementioned method by Ehinger and Lailly 
(1995)) by carrying out the demigration stage by means of ray tracing allowing 
calculation of multi-valuated arrivals. Such ray tracing is for example described in: 

Jurado, F., Lailly, P., and Ehinger, A., (1998), Fast 3D Two Point Raytracing 
5 forTravel Time Tomography: Proceedings of SPIE, Mathematical Methods in 
Geophysical Imaging V, 3453, 70-81. 

Then, after forming the various interface/velocity pairs by kinematic inversion 
at fixed velocity on the obtained time/ray parameters data, the kinematic information 
is calculated as described in the most general case by means, for example, of multi- 
10 offset ray tracing in each one of these models and for each multi-offset collection 
selected, but this time adapted to the multi-valuated arrival calculations as described 
in the aforementioned document by Jurado et al. (1998). Both prestack multi- 
valuated travel times and their associated ray parameters are thus obtained for each 
multi-offset collection selected. 

15 Then, the coherence analysis between the travel time curves thus predicted 

and the events considered in the seismic data is carried out as described in the most 
general case, but independently for each branch of the multi-valuated arrivals. The 
travel times thus collected for each branch of the multi-valuated arrivals are then 
injected, with their associated ray parameters, in a prestack kinematic inversion 

20 method allowing taking account of the multi-valuated arrivals, as described in the 
aforementioned document by Delprat-Jannaud and Lailly (1995). 

An application of the method to determination of a velocity model of seismic 
waves in an underground formation has been described. It is however obvious that 



the method can also apply to determination of the velocity of propagation of other 
types of waves in a heterogeneous model exhibiting discontinuities. 
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